Thresholding based segmentation¶

Motivation

Images can always contain a lot of information

Main character, objects, backgrounds, etc.

image info

To help computer identify different elements in an image, we use segmentation to divide image into sub-regions.

Definition

Thresholding based segmentation directly divide images into regions based on intensity values and/or properties of these values.

Suppose that the intensity histogram in this plot corresponds to an image

image info

To extract objects from background, one obvious way is to set a threshold that separates these modes.

Let f(x, y) being an image, and T being the threshold

If f(x, y) > T, these points are called object points

Otherwise, those points are called background points

In a segmented image, $g(x, y)= \begin{cases}1 & \text { if } f(x, y)>T \\ 0 & \text { if } f(x, y) \leq T\end{cases}$

There could also be some more difficult problems

image.png

In this case, the segmented image is given by $g(x, y)= \begin{cases}a & \text { if } f(x, y)>T_2 \\ b & \text { if } T_1 < f(x, y) \leq T_2 \\ c & \text { if } f(x, y) \leq T_1\end{cases}$

where a, b, and c are any three distinct intensity values.

Different Elements Influencing Thresholding¶

Having noise in the image can influence the histogram of an image image.png

Illumination and reflectance can also affect the histogram of an image and thus the quality of thresholding image.png

How To Find a Threshold¶

Basic Global Thresholding:

  1. Select an initial estimate for the global threshold, T.
  2. Segment the image using T. This will produce two groups of pixels: $G_1$ consisting of all pixels with intensity values > T, and $G_2$ consisting of pixels with values $\leq$ T.
  3. Compute the average (mean) intensity values $m_1$ and $m_2$ for the pixels in $G_1$ and $G_2$, respectively.
  4. Compute a new threshold value: $T=\frac{1}{2}(m_1+m_2)$
  5. Repeat step 1 to 4 until the difference between values of $T$ in successive iterations is smaller than a predefined parameter $\Delta T$

Finding Optimum: Limitations and solution

Thresholding can be viewed as a statistical-decision theory problem.

Objective: Minimize the average error incurred in assigning pixels to two or more groups.

Elegant closed-form solution: Bayes decision rule

Solution is based on:

  1. The probability density function (PDF) of the intensity levels
  2. Probability that each class occurs.

Issue: Estimating PDF is not a trivial task.

Alternative: Otsu's Method (Otsu [1979])

Why is it optimum?

Otsu's method maximizes the between-class variance, while variance is widely used in discriminant analysis.

Basic idea: Well-thresholded classes should be distinct with respect to the intensity values.

Therefore, such a threshold gives the best separation between classes.

Not only optimal, but also simple since entirely based on computations performed on this histogram.

The method can be summarized as follows:

  1. Compute the normalized histogram of the input image. Denote the components of the histogram by $p_i, i = 0, 1, 2,..., L - 1$
  2. Compute the cumulative sums, $P_{1}(k)$, for $k = 0, 1, 2,..., L - 1$.
  3. Compute the cumulative means, $m(k)$, for $k = 0, 1, 2, Á , L - 1$.
  4. Compute the global intensity mean $m_G$.
  5. Compute the between-class variance $\sigma_{B}^2(k)$ for $k = 0, 1, 2, Á , L - 1$.
  6. Obtain the Otsu threshold, $k*$, as the value of $k$ for which $\sigma^2_B(k)$ is maximum. If the maximum is not unique, obtain $k*$ by averaging the values of $k$ corresponding to the various maxima detected
  7. Obtain the separability measure, $\eta^*$.
In [33]:
import numpy as np
import os
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import skimage as ski
from skimage.color import rgb2gray
from PIL import Image
import cv2
import random
In [34]:
file_name = os.path.join('griez.jpeg')
our_image = ski.io.imread(file_name)
img = mpimg.imread(file_name)
imgplot = plt.imshow(img)
In [35]:
gs_image = ski.color.rgb2gray(img)
plt.imshow(gs_image, cmap = 'gray')
Out[35]:
<matplotlib.image.AxesImage at 0x7fd02ac13340>
In [36]:
n, bins, patches = plt.hist(gs_image.flatten(), bins=100)
In [37]:
var_list = []
for k in bins[1:]:
  p1 = 0
  p2 = 0
  for i in range(len(bins)-1):
    if bins[i] < k:
      p1 += n[i]/(np.shape(img)[0]*np.shape(img)[1])
    else:
      p2 += n[i]/(np.shape(img)[0]*np.shape(img)[1])
  mk = 0
  mk2 = 0
  for i in range(len(bins)-1):
    if bins[i] < k:
      mk += bins[i]*n[i]/(np.shape(img)[0]*np.shape(img)[1])
    else:
      mk2 += bins[i]*n[i]/(np.shape(img)[0]*np.shape(img)[1])
  var = p1*p2*(mk-mk2)**2
  var_list.append(var)

ix = var_list.index(max(var_list))
print(bins[ix])

0.38
In [38]:
for i in range(len(gs_image[:,0])):
  for j in range(len(gs_image[0,:])):
    if gs_image[i][j] < bins[ix]:
      gs_image[i][j] = 0
    else:
      gs_image[i][j] = 1
In [39]:
plt.imshow(gs_image, cmap = 'gray')
Out[39]:
<matplotlib.image.AxesImage at 0x7fcfc86379d0>

Multiple Thresholds¶

Thresholding techniques discussed before can be extended to arbitrary number of classes.

In the case of $K$ classes, $C_1, C_2, ..., C_K$m the between class variance can be expressed as

$\sigma_B^2=\sum_{k=1}^K P_k(m_k-m_G)^2$

For example, if we have two threshold, the between class variance is

$\sigma_B^2=P_1(m_1-m_G)^2+P_2(m_2-m_G)^2+P_3(m_3-m_G)^2$

Iteratively performing this method over the entire image will give us two optimal thresholds

In [40]:
file_name = os.path.join('griez.jpeg')
our_image = ski.io.imread(file_name)
img = mpimg.imread(file_name)
gs_image = ski.color.rgb2gray(img)
In [41]:
mg = 0
for i in range(len(bins)):
  if i < 100:
    mg += bins[i]*n[i]/(np.shape(img)[0]*np.shape(img)[1])

var_array = np.zeros((len(bins[1:]), len(bins[1:])))
var_list2 = []
max_idx = []
for k in range(len(bins[1:])):
  for k2 in range(len(bins[k+1:])):
    p1 = 0
    p2 = 0
    p3 = 0
    for i in range(len(bins)-1):
      if i <= k:
        p1 += n[i]/(np.shape(img)[0]*np.shape(img)[1])
      elif i <= k2:
        p2 += n[i]/(np.shape(img)[0]*np.shape(img)[1])
      else:
        p3 += n[i]/(np.shape(img)[0]*np.shape(img)[1])
    mk1 = 0
    mk2 = 0
    mk3 = 0
    for i in range(len(bins)-1):
      if i <= k:
        mk1 += bins[i]*n[i]/(np.shape(img)[0]*np.shape(img)[1])
      elif i <= k2:
        mk2 += bins[i]*n[i]/(np.shape(img)[0]*np.shape(img)[1])
      else:
        mk3 += bins[i]*n[i]/(np.shape(img)[0]*np.shape(img)[1])
    var = p1 * (mk1-mg)**2 + p2 * (mk2-mg)**2 + p3 * (mk3-mg)**2
    var_array[k,k2] = var
In [42]:
for i in range(len(var_array[0])):
  max_idx.append(np.max(var_array[i]))
In [43]:
np.where(var_array == max(max_idx))
Out[43]:
(array([38]), array([61]))
In [44]:
k1 = bins[38]
k2 = bins[61]
for i in range(len(gs_image[:,0])):
  for j in range(len(gs_image[0,:])):
    if gs_image[i][j] <= k1:
      gs_image[i][j] = 0
    elif k1 < gs_image[i][j] <= k2:
      gs_image[i][j] = 0.5
    else:
      gs_image[i][j] = 1
plt.imshow(gs_image, cmap ='gray')
Out[44]:
<matplotlib.image.AxesImage at 0x7fd0099f1700>

Variable Thresholding¶

One of the simplest approach: Image partition

Example:

  1. Original image
  2. Its histogram
  3. Using method 1 iteratively to obtain a threshold
  4. Using Otsu's method to obtain a threshold
  5. Partition
  6. Using Otsu's method for each sub-region and segmenting image-2.png
In [45]:
file_name = os.path.join('griez.jpeg')
our_image = ski.io.imread(file_name)
img = mpimg.imread(file_name)
gs_image = ski.color.rgb2gray(img)
In [46]:
from PIL import Image

def partition_image(image_path, num_rows, num_cols):
    
    image = Image.open(image_path)

    # Get the dimensions of the image
    width, height = image.size
    
    # Calculate the width and height of each rectangular region
    region_width = width // num_cols
    region_height = height // num_rows
    
    regions = []
    
    # Iterate over rows and columns to extract each region
    for row in range(num_rows):
        for col in range(num_cols):
            # Define the bounding box for the current region
            left = col * region_width
            upper = row * region_height
            right = left + region_width
            lower = upper + region_height
            
            # Crop the region from the image
            region = image.crop((left, upper, right, lower))
            
            # Append the region to the list of regions
            regions.append(region)
    
    return regions

# Example usage

image_path = './griez.jpeg'
num_rows = 2
num_cols = 2
regions = partition_image(image_path, num_rows, num_cols)
In [47]:
gs_region0 = ski.color.rgb2gray(regions[0])
gs_region1 = ski.color.rgb2gray(regions[1])
gs_region2 = ski.color.rgb2gray(regions[2])
gs_region3 = ski.color.rgb2gray(regions[3])
In [48]:
n_0, bins_0, patches_0 = plt.hist(gs_region0.flatten(), bins=100)
n_1, bins_1, patches_1 = plt.hist(gs_region1.flatten(), bins=100)
n_2, bins_2, patches_2 = plt.hist(gs_region2.flatten(), bins=100)
n_3, bins_3, patches_3 = plt.hist(gs_region3.flatten(), bins=100)
In [49]:
def otsu(n, bin):
  var_list = []
  for k in bin[1:]:
    p1 = 0
    p2 = 0
    for i in range(len(bin)-1):
      if bin[i] < k:
        p1 += n[i]/(np.shape(img)[0]*np.shape(img)[1])
      else:
        p2 += n[i]/(np.shape(img)[0]*np.shape(img)[1])
    mk = 0
    mk2 = 0
    for i in range(len(bin)-1):
      if bin[i] < k:
        mk += bin[i]*n[i]/(np.shape(img)[0]*np.shape(img)[1])
      else:
        mk2 += bin[i]*n[i]/(np.shape(img)[0]*np.shape(img)[1])
    var = p1*p2*(mk-mk2)**2
    var_list.append(var)

  ix = var_list.index(max(var_list))
  return bins[ix]
In [50]:
print(otsu(n_0, bins_0))
print(otsu(n_1, bins_1))
print(otsu(n_2, bins_2))
print(otsu(n_3, bins_3))

0.3
0.3
0.52
0.56
In [51]:
def plot(img, thres):
  for i in range(len(img[:,0])):
    for j in range(len(img[0,:])):
      if img[i][j] < thres:
        img[i][j] = 0
      else:
        img[i][j] = 1
  return img
In [52]:
img0 = plot(gs_region0, otsu(n_0, bins_0))
plt.imshow(img0, cmap = 'gray')
Out[52]:
<matplotlib.image.AxesImage at 0x7fd00b39a850>
In [53]:
img1 = plot(gs_region1, otsu(n_1, bins_1))
plt.imshow(img1, cmap = 'gray')
Out[53]:
<matplotlib.image.AxesImage at 0x7fcfc90a2c40>
In [54]:
img2 = ski.img_as_float32 = plot(gs_region2, otsu(n_2, bins_2))
plt.imshow(img2, cmap = 'gray')
Out[54]:
<matplotlib.image.AxesImage at 0x7fcfd06b18b0>
In [55]:
img3 = plot(gs_region3, otsu(n_3, bins_3))
plt.imshow(img3, cmap = 'gray')
Out[55]:
<matplotlib.image.AxesImage at 0x7fd00af59880>
In [56]:
width, height = np.shape(img0)
empty_image = np.zeros((width*2, height*2))
for i in range(0, width*2):
  for j in range(0, height*2):
    if i < 333 and j >= 500:
      empty_image[i,j] = img1[i,j-500]
    if i < 333 and j < 500:
      empty_image[i,j] = img0[i,j]
    if i >= 333 and j >= 500:
      empty_image[i,j] = img3[i-333,j-500]
    if i >= 333 and j < 500:
      empty_image[i,j] = img2[i-333,j]
In [57]:
plt.imshow(empty_image, cmap = 'gray')
Out[57]:
<matplotlib.image.AxesImage at 0x7fcff9443640>

Application¶

White blood cell segmentation from microscopic blood images¶

Methods:

  1. Pre-processing: Median filter are applied to the image to remove noise, and histogram equalization is applied on the green channel of the RGB input image to enhance the image quality.
  2. Segmentation: Goal is to segment the white blood cells (WBCs) from the other components of the image (i.e. red blood cells (RBCs) and background). Thresholding segmentation is applied. Both global and local thresholding are considered here.
  3. Image cleaning: Resulted image contains unwanted objects. Morphological opening operator is applied.

image.png

image.png

image.png

Reference: [1]. Gonzalez, Rafael C., and Richard E. Woods. Digital Image Processing. Pearson, 2019. [2]. Mohammed, Zhana Fidakar, and Alan Anwer Abdulla. “Thresholding-based white blood cells segmentation from microscopic blood images.” UHD Journal of Science and Technology, vol. 4, no. 1, 13 Feb. 2020, pp. 9–17, https://doi.org/10.21928/uhdjst.v4n1y2020.pp9-17.